﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Text.RegularExpressions;
using System.ComponentModel.Composition;

namespace Nature.Numerics
{
    /// <summary>
    /// Represents an algorithm for solving systems of linear algebraic equations
    /// </summary>
    public interface ILinearAlgebraicEquationsSolver
    {
        /// <summary>
        /// Solves the given system of linear algebraic equations. A return value indicates whether the algorithm succeeded.
        /// </summary>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <returns>true if the system was solved successfully; otherwise, false.</returns>        
        bool TrySolve(double[][] a, double[] b);
    }

    /// <summary>
    /// The exception that is thrown when solving a system of linear algebraic equations has failed
    /// </summary>
    [Serializable]
    public class FailedSolvingLinearAlgebraicEquations : System.Exception 
    {
        /// <summary>
        /// Initializes a new instance of the FailedSolvingLinearAlgebraicEquations class.
        /// </summary>
        public FailedSolvingLinearAlgebraicEquations() { }

        /// <summary>
        /// Initializes a new instance of the FailedSolvingLinearAlgebraicEquations class with a specified error message.
        /// </summary>
        /// <param name="message">The error message that explains the reason for the exception.</param>
        public FailedSolvingLinearAlgebraicEquations(string message) { }

        /// <summary>
        /// Initializes a new instance of the FailedSolvingLinearAlgebraicEquations class with a specified error message 
        /// and a reference to the inner exception that is the cause of this exception.
        /// </summary>
        /// <param name="message">The error message that explains the reason for the exception.</param>
        /// <param name="exception">The exception that is the cause of the current exception. 
        /// If the innerException parameter is not a null reference (Nothing in Visual Basic), 
        /// the current exception is raised in a catch block that handles the inner exception.</param>
        public FailedSolvingLinearAlgebraicEquations(string message, Exception exception) { }


    }

    public static partial class Extensions
    {
        /// <summary>
        /// Solves the given system of linear algebraic equations.
        /// </summary>
        /// <typeparam name="TLinearAlgebraicEquationsSolver">Genetic type implementing the cref="Nature.Numerics." ILinearAlgebraicEquationsSolver interface</typeparam>
        /// <param name="solver">Instance of the genetic type implementing the ILinearAlgebraicEquationsSolver interface</param>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <exception cref="Nature.Numerics.FailedSolvingLinearAlgebraicEquations"></exception>
        public static void Solve<TLinearAlgebraicEquationsSolver>(this TLinearAlgebraicEquationsSolver solver,
            double[][] a, double[] b)
            where TLinearAlgebraicEquationsSolver : ILinearAlgebraicEquationsSolver
        {
            if (false == solver.TrySolve(a, b))
            {
                throw new FailedSolvingLinearAlgebraicEquations();
            }
        }

        /// <summary>
        /// Solves the given system of linear algebraic equations. If the system was solved successfully, 
        /// returns an array of values representing the solution vector (x-vector); otherwise, null.
        /// </summary>
        /// <typeparam name="TLinearAlgebraicEquationsSolver">Genetic type implementing the ILinearAlgebraicEquationsSolver interface</typeparam>
        /// <param name="solver">Instance of the genetic type implementing the ILinearAlgebraicEquationsSolver interface</param>
        /// <param name="aString">A string representation of the system coefficients (a-matrix)</param>
        /// <param name="bString">A string representation of the constant terms (b-vector)</param>        
        /// <returns>An array of values representing the solution vector (x-vector), if the system was solved successfully, otherwise, null</returns>
        public static double[] TrySolve<TLinearAlgebraicEquationsSolver>(this TLinearAlgebraicEquationsSolver solver,
            string aString, string bString)
            where TLinearAlgebraicEquationsSolver : ILinearAlgebraicEquationsSolver
        {
            double[][] matrix = NumericalTypeConverter.ConvertToMatrix(aString);
            double[] array = NumericalTypeConverter.ConvertToArray(bString);
            if (matrix.Length != array.Length)
            {
                throw new ArgumentException();
            }
            if (solver.TrySolve(matrix, array))
            {
                return array;
            }
            return null;
        }

        /// <summary>
        /// Solves the given system of linear algebraic equations.
        /// </summary>
        /// <typeparam name="TLinearAlgebraicEquationsSolver">Genetic type implementing the ILinearAlgebraicEquationsSolver interface</typeparam>
        /// <param name="solver">Instance of the genetic type implementing the ILinearAlgebraicEquationsSolver interface</param>
        /// <param name="aString">A string representation of the system coefficients (a-matrix)</param>
        /// <param name="bString">A string representation of the constant terms (b-vector)</param>
        /// <returns>An array of values representing the solution vector (x-vector)</returns>
        /// <exception cref="Nature.Numerics.FailedSolvingLinearAlgebraicEquations"></exception>
        public static double[] Solve<TLinearAlgebraicEquationsSolver>(this TLinearAlgebraicEquationsSolver solver,
            string aString, string bString)
            where TLinearAlgebraicEquationsSolver : ILinearAlgebraicEquationsSolver
        {
            double[][] matrix = NumericalTypeConverter.ConvertToMatrix(aString);
            double[] array = NumericalTypeConverter.ConvertToArray(bString);
            if (matrix.Length != array.Length)
            {
                throw new ArgumentException();
            }
            if (solver.TrySolve(matrix, array)) { return array;  }
            throw new FailedSolvingLinearAlgebraicEquations();
        }
    }

    /// <summary>
    /// AHRpGauss2004 realization of the Gaussian Elimination algorithm for solving systems of linear algebraic equations.
    /// </summary>
    [Export(typeof(ILinearAlgebraicEquationsSolver))]
    public class AHRpGauss2004 : ILinearAlgebraicEquationsSolver
    {
        #region Static Members

        /// <summary>
        /// Solves the given system of linear algebraic equations. A return value indicates whether the algorithm succeeded.
        /// </summary>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <returns>true if the system was solved successfully; otherwise, false.</returns>    
        public static bool TrySolve(double[][] a, double[] b)
        {
            int nva = b.Length;
            for (int i = 0; i < nva - 1; ++i)
            {
                #region Pivoting
                {
                    int imax = i;
                    double scale = 0.0;
                    for (int ii = i; ii < nva; ++ii)
                    {
                        double[] rr = a[ii];
                        double current, maximum = 0.0;
                        for (int j = i; j < nva; ++j)
                        {
                            current = rr[j];
                            if (current < 0.0) { current = (-current); }
                            if (maximum < current) { maximum = current; }
                        }
                        if (maximum == 0.0)
                        {
                            return false;
                        }
                        current = rr[i] / maximum;
                        if (current < 0.0) { current = (-current); }
                        if (scale < current)
                        {
                            scale = current;
                            imax = ii;
                        }
                    }
                    if (imax != i)
                    {
                        double[] atmp = a[i];
                        a[i] = a[imax];
                        a[imax] = atmp;

                        double btmp = b[i];
                        b[i] = b[imax];
                        b[imax] = btmp;
                    }
                }
                #endregion

                double[] mainRow = a[i];
                double mainElement = mainRow[i];
                for (int ii = i + 1; ii < nva; ++ii)
                {
                    double[] row = a[ii];
                    double factor = row[i] / mainElement;
                    b[ii] -= b[i] * factor;
                    for (int j = i + 1; j < nva; ++j)
                    {
                        row[j] -= mainRow[j] * factor;
                    }
                }
            }


            for (int i = nva - 1; i >= 0; --i)
            {
                double[] row = a[i];
                for (int j = i + 1; j < nva; j++)
                {
                    b[i] -= b[j] * row[j];
                }
                b[i] /= row[i];
            }
            return true;
        }

        /// <summary>
        /// Solves the given system of linear algebraic equations.
        /// </summary>
        /// <param name="solver">Instance of the genetic type implementing the ILinearAlgebraicEquationsSolver interface</param>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <exception cref="Nature.Numerics.FailedSolvingLinearAlgebraicEquations"></exception>
        public static void Solve(double[][] a, double[] b)
        {
            if (false == AHRpGauss2004.TrySolve(a, b))
            {
                throw new FailedSolvingLinearAlgebraicEquations();
            }
        }

        #endregion

        #region ILinearAlgebraicEquationsSolver Members

        /// <summary>
        /// Solves the given system of linear algebraic equations. A return value indicates whether the algorithm succeeded.
        /// </summary>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <returns>true if the system was solved successfully; otherwise, false.</returns>  
        bool ILinearAlgebraicEquationsSolver.TrySolve(double[][] a, double[] b) { return AHRpGauss2004.TrySolve(a, b); }

        #endregion        
    }

    /// <summary>
    /// AHCpGauss2004 realization of the Gaussian Elimination algorithm for solving systems of linear algebraic equations.
    /// </summary>
    [Export(typeof(ILinearAlgebraicEquationsSolver))]
    public class AHCpGauss2004 : ILinearAlgebraicEquationsSolver
    {
        #region Nested Types

        /// <summary>
        /// A set of data necessary to carry out the computations
        /// </summary>
        struct Data
        {
            double[][] m_a;
            double[] m_b;
            int[] m_columns;
            int m_nva;

            public Data(int nva)
                : this()
            {
                m_columns = new int[nva];
                m_b = new double[nva];
                m_nva = nva;
            }

            public void Import(double[][] a, double[] b)
            {
                m_a = a;
                for (int j = 0; j < m_b.Length; ++j)
                {
                    m_columns[j] = j;
                    m_b[j] = b[j];
                }
            }

            public void Export(double[] b)
            {
                for (int j = 0; j < m_nva; ++j)
                {
                    int index = m_columns[j];
                    b[index] = m_b[j];
                }
            }

            public double this[int i]
            {
                get { return m_b[i]; }
                set { m_b[i] = value; }
            }

            public double this[int i, int j]
            {
                get { return m_a[i][m_columns[j]]; }
                set { m_a[i][m_columns[j]] = value; }
            }

            public void SwapColumns(int j1, int j2)
            {
                int jr = m_columns[j1];
                m_columns[j1] = m_columns[j2];
                m_columns[j2] = jr;
            }
        }

        /// <summary>
        /// Encapsulates a set of data necessary to carry out the computations
        /// </summary>
        class Provision
        {
            /// <summary>
            /// A set of data necessary to carry out the computations
            /// </summary>
            public Data Data { get; set; }
        }

        #endregion

        /// <summary>
        /// Solves the given system of linear algebraic equations. A return value indicates whether the algorithm succeeded.
        /// </summary>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <param name="provisionObject">An object encapsulating a set of data necessary to carry out the computations</param>
        /// <returns>true if the system was solved successfully; otherwise, false.</returns>
        public bool TrySolve(double[][] a, double[] b, ref object provisionObject)
        {
            int nva = a.Length;
            Provision provision = provisionObject as Provision;
            if (provision == null)
            {
                provision = new Provision { Data = new Data(nva) };
                provisionObject = provision;
            }
            Data data = provision.Data;
            data.Import(a, b);
            for (int i = 0; i < nva - 1; ++i)
            {
                #region Pivoting
                {
                    double current, maximum = 0.0;
                    int jmax = i;
                    for (int j = i; j < nva; ++j)
                    {
                        current = data[i, j];
                        if (current < 0.0) { current = (-current); }
                        if (maximum < current)
                        {
                            jmax = j;
                            maximum = current;
                        }
                    }
                    if (maximum == 0.0) { return false; }
                    data.SwapColumns(i, jmax);
                }
                #endregion

                double mainElement = data[i, i];
                for (int ii = i + 1; ii < nva; ++ii)
                {
                    double factor = data[ii, i] / mainElement;
                    data[ii] -= data[i] * factor;
                    for (int j = i + 1; j < nva; ++j)
                    {
                        data[ii, j] -= data[i, j] * factor;
                    }
                }
            }


            for (int i = nva - 1; i >= 0; --i)
            {
                for (int j = i + 1; j < nva; j++)
                {
                    data[i] -= data[j] * data[i, j];
                }
                data[i] /= data[i, i];
            }

            data.Export(b);
            return true;
        }

        #region ILinearAlgebraicEquationsSolver Members

        /// <summary>
        /// Solves the given system of linear algebraic equations. A return value indicates whether the algorithm succeeded.
        /// </summary>
        /// <param name="a">Coefficients of the system (a-matrix)</param>
        /// <param name="b">Constant terms (b-vector)</param>
        /// <returns>true if the system was solved successfully; otherwise, false.</returns>
        public bool TrySolve(double[][] a, double[] b)
        {
            object provision = null;
            return this.TrySolve(a, b, ref provision);
        }

        #endregion
    }
}
